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Abstract 

We propose an algorithm to estimate the Hurst exponent of high-dimensional fractals, based 
on a generalized high-dimensional variance around a moving average low-pass filter. As working 
examples, we consider rough surfaces generated by the Random Midpoint Displacement and by the 
Cholesky-Levinson Factorization algorithms. The surrogate surfaces have Hurst exponents ranging 
from 0.1 to 0.9 with step 0.1, and different sizes. The computational efficiency and the accuracy 
of the algorithm are also discussed. 

PACS numbers: 05.10.-a, 05.40.-a, 05.45.Df, 68.35.Ct 



I. INTRODUCTION 



The scahng properties of random curves and surfaces can be quantified in terms of the 
Hurst exponent iJ, a parameter defined in the framework of the fractional Brownian walks 
introduced in [1]. A fractional Brownian function f{r) : M'^ ^ M, is characterized by a 
variance aj^: 

al= ([/(r + A)-/(r)]2) oc ||A|r with a = 2H , (1) 

with r = {xi, X2, Xd) , A = (Ai, A2, A^) and ||A|| = ^y Xf + + ... + X'^ ; a power 
spectrum Sh- 

Sh oc \\u:\\-'^ with (3 = d + 2H , (2) 

with ci; = {uJi,uj2, ■■■,^d) the angular frequency, \\(jl>\\ = a/o;^ + a;| + . . . + tu^ ; a number of 
objects Nh of characteristic size e needed to cover the fractal: 

Nh oc e"^ with D = d+1 - H , (3) 

D being the fractal dimension of f{r). The Hurst exponent ranges from to 1, taking 
the values H = 0.5, H > 0.5 and H < 0.5 respectively for uncorrelated, correlated and 
anticorrelated Brownian functions. 

The application of fractal concepts, through the estimate of H, has been proven useful 
in a variety of fields. For example in d = 1, heartbeat intervals of healthy and sick hearts 
are discriminated on the basis of the value of if [21 E] ; the stage of financial market devel- 
opment is related to the correlation degree of return and volatility series |1] ; coding and non 
coding regions of genomic sequences have different correlation degrees [5]; climate models 
are validated by analyzing long-term correlation in atmospheric and oceanographic series 
P [7] . In > 2 fractal measures are used to model and quantify stress induced morpholog- 
ical transformation ^J; isotropic and anisotropic fracture surfaces P, [THl El [111 HSj; static 
friction between materials dominated by hard core interactions ; diffusion [16] and 
transport [TTIIIH] in porous and composite materials; mass fractal features in wet/dried gels 
[I9] and in physiological organs (e.g. lung) [20]; hydrophobicity of surfaces with hierarchic 
structure undergoing natural selection mechanism and solubility of nanoparticles |22]; 
digital elevation models [23] and slope fits of planetary surfaces [2^ . 
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A number of fractal quantification methods - based on the Eqs. ([lj|3j) or on variants of 
these relationships - like Rescaled Range Analysis (R/S), Detrended Fluctuation Analysis 
(DFA), Detrending Moving Average Analysis (DMA), Spectral Analysis, have been thus 
proposed to accomplish accurate and fast estimates of H in order to investigate correlations 
at different scales m d = 1. A comparatively small number of methods able to capture spa- 
tial correlations-operating m. d> 2-have been proposed so far [251 ESI ETJ EU [29l [SQl |3T1 132] . 
This work is addressed to develop an algorithm to estimate the Hurst exponent of high- 
dimensional fractals and thus is intended to capture scaling and correlation properties over 
space. The proposed method is based on a generalized high-dimensional variance of the 
fractional Brownian function around a moving average. In Section [IT], we report the rela- 
tionships holding for fractals with arbitrary dimension. It is argued that the implementation 
can be carried out in directed or isotropic mode. We show that the Detrending Moving Av- 



erage (DMA) method [30l |3T| [32] is recovered for = 1. In Section III, the feasibility of 
the technique is proven by implementing the algorithm on rough surfaces - with different 
size A^^i X N2 and Hurst exponent H - generated by the Random Midpoint Displacement 
(RMD) and by the Cholesky-Levinson Factorization (CLF) methods [33l IM] . The gener- 
alized variance is estimated over sub-arrays ni x n2 with different size ( "scales") and then 
averaged over the whole fractal domain A^i x N2. This feature reduces the bias effects due 
to nonstationarity with an overall increase of accuracy - compared to the two-point corre- 
lation function, whose average is calculated over all the fractal. Furthermore - compared to 
the two-point correlation function, whose implementation is carried out along 1-dimensional 
lines (e.g. for the fracture problem, the two-point correlation functions are measured along 
the crack propagation direction and the perpendicular one), the present technique is carried 



out over d-dimensional structures (e.g. squares va. d = 2). In Section IV, we discuss accuracy 
and range of applicability of the method. 

II. METHOD 

In order to implement the algorithm, the generalized variance uI^ma introduced: 



Afi~mi N2-m2 Na-ma 

M . 



^'dMA = XI 5Z 5Z [/(^1.^2,-.,^d) - 7ni,n2,...n,(«l,«2,...,id) , (4) 
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where f{ii,i2,...,id) = f{i) is a fractional Brownian function defined over a discrete d- 
dimensional domain, with maximum sizes Ni, N2, N^. It is ii = 1,2,..., Ni, 12 = 
1, 2, N2, id = 1, 2, Nd. n = {ni,n2, Ud) defines the sub-arrays z/^ of the fractal 
domain with maximum values nimax = max{ni}, n2max = max{n2}, Udmax = niax{nrf}; 
mi = mt{ni9i) , 1712 = int{n292),...,md = mt{nd9d) and 9i, 92, ...9d are parameters 
ranging from to 1; A/" = (A^i - nimax) ■ {N2 - n2max) ■ ■■■ ■ {Nd - Udmax)- The function 
7ni,n2,...,nd(H,^2, = 7 is giveu by: 

^ ni — l—nii 112 — I— m2 

/ni,n2,...,nti(^l' ^2) ^d) — / ^ / ^ 

fci=— mi k2=—m2 

Ud-l-ma 

... ^ f{ii-ki,i2-k2,...,id-kd), (5) 

fcd=-m,d 

that is an average of / calculated over the sub-arrays Vd- The Eqs. Q and Q are defined 
for any value of ni,n2, Ud and for any shape of the sub-arrays, however, it is preferable 
to choose sub-arrays with rii = n2 = .... = rid to avoid spurious directionality in the results. 
The generahzed variance cr|,j^^^ varies as + + ... + )'^^ as a consequence of the 

property ([T| of the fractional Brownian functions. 

Upon variation of the parameters 9i, 92 , ... ,9d in the range [0, 1], the indexes 11,12, id 
and ki, k2,...,kd of the sums in the Eqs. Q and ([S]) are accordingly set within Ud- In 
particular, {ii, 12, id) coincides respectively with: (a) one of the vertices of Vd for 9i = 92 = 
... = 9d = and 9i = 92 = ... = 9d = 1 or (b) the center of Ud for 9i = 92 = ... = 9d = 1/2. 
It is worthy of note that the choice 6*1 = ^2 = ••• = dd = 1/2 corresponds to the isotropic 
implementation of the algorithm, while 9i = 92 = ... = 9d = ^ and 9i = 92 = ... = 
9d = 1 correspond to the directed implementation. For example in d = 2, the isotropic 
implementation implies that the variance defined by the Eq. Q is referred to a moving 
average / calculated over squares ni x n2 whose center is {ii,i2)- Conversely, the directed 
implementation implies that the function / is calculated over squares ni x n2 with one of the 
vertices in {ii, i2). The directed mode is of interest to estimate H in fractals with preferential 
growth direction, e.g. biological tissues (lung), epitaxial layers, crack propagation in fracture 
[anisotropic fractals). If the fractal is isotropic and the accuracy is a priority, the parameters 
9i, 92, 9d should be preferably taken equal to 1/2 to achieve the most precise estimate of 
H. The dependence of the algorithm on 6* for d = 1 has been discussed in [32] . 



In order to calculate the Hurst exponent, the algorithm is implemented through the 
following steps. The moving average / is calculated for different sub-arrays z/^, by 
varying ni,n2,..., from 2 to the maximum values nimax,n2max, ridmax- The val- 
ues nimax,n2max, ^dmax depend on the maximum size of the fractal domain. In or- 
der to minimize the saturation effects due to finite-size, it should be: riimax << Ni] 



n2max « ^2'-, Udmax « Nd- Thesc coustraiuts will be further clarified in Section [ITT 
where the algorithm is implemented over fractal surfaces with different sizes. For each 
sub- array Ud, the corresponding value of Cdma calculated and finally plotted on log-log 
axes. 

To elucidate the way the algorithm works, in the following we consider its implementation 
for d = 1 and d = 2. The case d = 1 reduces to the Detrending Moving Average (DMA) 
method already used for long-range correlated time series [301 EH |32] • 

1- dimensional case: 

By posing c? = 1 in the Eq. Q, one obtains: 



2 

^DMA 



Ni-nii ^ 2 



Nl — riimax 



E ' (6) 



where Ni is the length of the sequence, rii is the sliding window and ni^ax = maxjni} ^ A^^i. 
The quantity uii = int(r;,i^^i) is the integer part of ni6i and 6i is a parameter ranging from 
to 1. The relationship ^ defines a generalized variance of the sequence f{ii) with respect 
to the function fmih)'- 



^1 , 

fel=— mi 

which is the moving average of f{ii) over each sliding window of length ui. The moving 
average /^(h) is calculated for different values of the window ni, ranging from 2 to the 
maximum value ni^ax- The variance (jj^MA then calculated according to the Eq. ^ and 
plotted as a function of ni on log-log axes. The plot is a straight line, as expected for a 
power-law dependence of on ni. 

^Ima ~ < ■ (8) 
The Eq. (fsl) allows one to estimate the scaling exponent H of the series f{ii). Upon variation 
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of the parameter 6i in the range [0, 1], the index ki in fmih) is accordingly set within the 
window Til. In particular, 9i = corresponds to average /„j (ii) over all the points to the left 
of ii within the window ni, 6i = 1 corresponds to average fmih) over all the points to the 
right of ii within the window rii, 6i = ^ corresponds to average fmih) with the reference 
point in the center of the window rii. 

2- dimensional case 

For (i = 2, the generalized variance defined by the Eq.Q writes: 



2 

'-'DMA 



{Ni — nimax){N2 — n2max) 



with fm,n2{ii,i2) given by: 



Ni-mi N2-m2 ^ 

/(^I'^a) - 7ni,n2(n,«2) , (9) 



«i=rii-mi i2=n2-m2 



ni — 1— mi n2 — 1— m2 



nin2 



Yl Y /(^i - ^i'^2 - ^2) 



(10) 



k\=—mi k2=—m2 

The average / is calculated over sub-arrays with different size ni x ^2. The next step is the 
calculation of the difference f{ii, 12) — fni,n2{H, ^2) for each sub-array rii x n2- A log-log plot 



DMA 



2H 



~ s 



H 



(11) 



as a function of s = nf + n"^, yields a straight line with slope H. 

Depending upon the values of the parameters 61 and 62, entering the quantities mi = 
int(ni6'i) and m2 = int(n26'2) in the Eqs. ( 9p0 ), the position of {ki,k2) and {ii,i2) can be 
varied within each sub-array. {ii,i2) coincides with a vertex of the sub-array if: (i) 61 = 0, 
62 = 0; {ii) 61 = 0, 6*2 = 1; {Hi) 61 = 1, 62 = 0; {iv) 61 = 1, 62 = 1 {directed implementation) . 
The choice 6'i = ^^2 = 1/2 corresponds to take the point {ii,i2) coinciding with the center of 
each sub-array rii x n2 {isotropic implementation) [41j . 



III. RESULTS 

In order to test feasibility and robustness of the proposed method, synthetic rough surfaces 
with assigned Hurst exponents have been generated by the Random Midpoint Displacement 
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(RMD) algorithm and by the Cholesky-Levinson Factorization (CLF) method [331 [31]. The 
widespread use of the RMD algorithm is based on the trade-off of its fast, simple and efficient 
implementation to its limited accuracy especially for H <^ 0.5 and H ^ 0.5. Conversely, the 
Cholesky-Levinson Factorization method is one of the most accurate techniques to generate 
Id and 2d fractional Brownian functions, at the expenses of a more complex implementation 
structure [42J. 

In Fig. [T| the log-log plots of crj^MA ^ function of s are shown for the synthetic fractal 
surfaces generated by the RMD (circles) and by the CLF method (squares). The surfaces 
have Hurst exponents iJj„ ranging from 0.1 to 0.9 with step 0.1. The domain sizes are 
respectively A^i x A^2 = 256 x 256 (a), A^i x A^2 = 1024 x 1024 (b) and NixN2 = 4096 x 4096 
(c). The dashed lines show the behavior that should be exhibited by variances varying 
exactly as s^'"^ over the entire range of scales. The plots of o"!,^^^ as a function s are in good 



agreement with the behavior expected on the basis of the Eq. ( 11 ). The quality of the fits is 
higher for the surfaces generated by the CLF method, confirming that the RMD algorithm 
synthesizes less accurate fractals. By comparing the results of the simulation (symbols) to 
the straight lines corresponding to full linearity over the whole range (dashed), deviations 
from the full linearity can be observed especially for the small surfaces at the extremes of 
the scale. A plot of the slopes for the fractal surfaces generated by the CLF algorithm is 
shown in Fig. [2] for different sizes of the fractal domain. A detailed discussion of the origin 



of the deviations at low and large scales is reported in the Section IV 



Finally, we also show three examples of digital images currently mapped to fractal surfaces 
with reference to the color intensity i.e. to the levels of Red, Green and Blue (RGB). The 
Hurst exponents estimated by the proposed method are respectively H = 0.1 (a), H = 0.5 
(b) and H = 0.9 (c) for the images in Fig. |3j 



IV. DISCUSSION 



The proposed algorithm is characterized by short execution time and ease of implemen- 
tation. By considering the case d = 2, the function /ni,n2(^i) "^2) is indeed simply obtained 
by summing the values of f{ii,i2) over each sub-array ni x ^2. Then the sum is updated 
at each step by adding the last and discarding the first row (column) of each sliding array 
111 X n2. For higher dimensions, the sum is updated at each step by adding and discarding a 
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FIG. 1: (Color online) Log-log plot of cr'jjj^^^ for fractal surfaces respectively with size A'^i x N2 = 
256 X 256 (a), iVi x iV2 = 1024 x 1024 (b) and iVi x iV2 = 4096 x 4096 (c). The data refer to fractal 
surfaces generated by the RMD (circles) and by the CLF (squares) methods. The Hurst exponent 
Hin - input of the RMD and the CLF algorithm - varies from 0.1 to 0.9 with step 0.1. The results 
correspond to the isotropic implementation, i.e. with the parameters = ^2 = 1/2 in the Eq.Q. 
The dashed lines represent the behavior expected for full linearity, i.e. the log-log plot of curves 
varying as s^"\ It is worthy of note that the CLF data are closer to the full-linearity compared 
to the RMD ones. 

d — 1 dimensional set of each array ni x 71,2 x ... x n^. The algorithm does not use arbitrary 
parameters, the computation simply relying on averages of /. We will now argue on the 
origin of the deviations at small and at large scales. 

Deviations at large scales. The deviations from the linearity at large scales, leading to 
the saturation of the are due to finite size effects. The small surfaces do not contain 

enough data to make the evaluation of the scaling law over the sub-arrays statistically 
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FIG. 2: (Color online) Plot of the values of H obtained by linear fit of the curves shown in 
Fig. [T] (a), (b), (c). The data refer to the fractal surfaces generated by the Cholesky-Levinson 
Factorization method (squares). The dashed lines represent the ideal behavior: H = Hm- 




FIG. 3: (Color online). Cloudy sky images respectively with Hurst exponent H = 0.1 (a), H = 0.5 
(b) and H = 0.9 (c). Such heterogeneous images are represented as fractal surfaces by mapping 
the color intensity (RGB content). 

meaningful. By comparing the data in Figs, [l] (a), (b), (c), one can note that the saturation 
effect decreases upon increasing the size A^^i x N2 of the fractal surface. The finite size effects 
become negligible when the conditions riimax « Ni, n2max « N2] ridmax « Nd are 
fulfilled. 

Deviations occurring at small scales. The deviations occurring at low scales are related 
to the departure of the low-pass filter from the ideality. This problem also occurs with one- 
dimensional fractals (time series) resulting in the quite generally reported overestimation of 
H in anticorrelated signals and underestimation of H in correlated signals [351 EHl EH EH |39] . 
We will discuss the origin of these deviations by means of the filter transfer function T-Cr{uj) 
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[40j . The algorithm is based on a generahzed variance of the function / with respect to 
/. The function / is the output of a low-pass filter driven by /, with impulse response a 
box-car function. In the Appendix, the transfer function Tiriuj) of / is explicitly calculated 
and shown in Fig.Q for d = 2. For an ideal low-pass filter, the transfer function should be 
one or zero respectively at frequencies lower or higher than the cut-off frequency. However, 
in real low-pass filters, at frequencies lower than the cut-off frequency, all the components 
of the signal suffer some attenuation but u = 0. The cut-off frequencies of Tirioj) are 



Ui = vr/rj, i.e. the first zeroes of the functions sincUjTj/cUjrj in the Eq. (A. 4). Moreover, in 
real filters, at frequencies higher than vr/xj, due to the presence of the sidelobes, components 
of the signals lying in the bands (vr/rj, 27r/rj); (27r/rj, Sn/ri); are not fully filtered out. As 
a result, the function / contains: (a) less components with frequency lower than Ui = n/ri 
and (b) more components with frequency higher than uJi = tt/ti compared to what it would 
be expected with an ideal low-pass filter. The lack of low-frequency components depends 
on the central lobe, while the excess of high-frequency components depends on the side 
lobes. The excess of high-frequency components results in a smaller value of the difference 
/ — /, i.e. in a decrease of (jj^MA ^T^d, thus, in an increase of the slope of the log-log plot. 
Conversely, the lack of low-frequency components results in a larger value of the difference 
/ — /, i.e. in an increase of crjjj^^ and, thus, in a decrease of the slope of the log-log plot. The 
two effects are more relevant with smaller values of the scales, when the filter nonideality 
is greater. Moreover, as one can deduce from the Eqs. ^ and (A.6), the effect of the side 
lobes dominates in high-frequency rich fractals with H < 0.5, while the effect of the central 
lobe is dominant in fractals with H > 0.5, rich of low-frequency components. 

In order to gain further insight in the above theoretical arguments, we report in Table 
|l]the slopes Hj, Hu and Hui of the curves (squares) plotted in Fig. [l] (b)) over different 
ranges. The slopes have been calculated by linear fit respectively over the ranges 10 < 
s < 100 {Hi), 10 < s < 1000 [Hn) and 10 < s < 10000 {Hm). The relative errors 
/S.H = (H — Hi„)/Hin are given respectively in the 3^'^, 5^^ and 7*'^ columns. The slope Hj 
is greater than the expected value Hm. The slope Hjj is overestimated for H = 0.1 and 
H = 0.2 and underestimated for H > 0.2. The slope Hjjj is underestimated since the effects 
of the finite-size of the fractal domain play a dominant role. 

We address the question if the artifacts due to the filter nonideality described above 
might be corrected somehow. In the remaining of this section, we will thus consider the use 
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FIG. 4: (Color online) Plot of the transfer function 7i'T{oJi,i^2)- 



TABLE I: Slopes Hm and relative errors AHj, AHu, AHjjj of the curves plotted 

in Fig. [l]^b) (squares). The slopes have been calculated by linear fit respectively over the ranges: 
10 < s < 100 (Hi), 10 < s < 1000 (Hn) and 10 < s < 10000 {Hm). The errors AHj, AHu, 



AHni 


are calculated 


as AH 


= {H- 


Hin ) / Hin • 


















Hi 


AHi 


Hn 


AHjj 




Hni 


AHni 




0.1 


0.1346 


+3.46 


10-1 


0.1073 


+7.30 


10" 


-2 


0.0718 


-2.822 


10" 


-1 


0.2 


0.2272 


+1.36 


10-1 


0.2050 


+2.50 


10" 


-2 


0.1700 


-1.500 


10- 


-1 


0.3 


0.3233 


+7.77 


10-2 


0.2995 


-1.67 


10" 


-3 


0.2716 


-9.467 


10- 


-2 


0.4 


0.4205 


+5.12 


10-2 


0.3970 


-7.50 


10" 


-3 


0.3691 


-7.725 


10- 


-2 


0.5 


0.5178 


+3.56 


10-2 


0.4973 


-5.40 


10" 


-3 


0.4752 


-4.960 


10- 


-2 


0.6 


0.6171 


+2.85 


10-2 


0.5973 


-4.50 


10- 


-3 


0.5617 


-6.383 


10- 


^2 


0.7 


0.7185 


+2.64 


10-2 


0.6956 


-6.29 


10- 


-3 


0.6770 


-3.286 


10- 


^2 


0.8 


0.8207 


+2.58 


10-2 


0.7999 


-1.25 


10" 


-4 


0.7659 


-4.263 


10- 


^2 


0.9 


0.9253 


+2.81 


10-2 


0.8999 


-1.11 


10" 


-4 


0.8679 


-3.567 


10" 


-2 



of windows whose general effect is to increase the width of the central lobe while reducing 
those of the sidelobes of the function Tiriuj) (a detailed description of these methods can 
be found in [®J). By restricting our discussion to the present technique, the correction is 
performed by using the following variant of the relationship (Is]) : 
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FIG. 5: (Color online) Plot of the function with / defined by the Eq. Ml) (solid lines) and 



/* defined by the Eq. ( 12 ) (dashed lines). It can be noted that the deviations of the slope at small 
scales are reduced by the use of /* implying a corresponding reduction of the relative error Ai7/ 
of two orders of magnitude. 



/ni,n2,...,n,(^1.^2,--, id) = (1 - «) /„i „^ (il , ^2 , • • • , Id) 

+a/ni,n2,...,nd(^l " 1 , ^2 " 1, • • • , «d " 1) 



(12) 



where a = nin2...nd/[{ni + l){n2 + l)...{nd + 1)]. The Eq. (12) reduces for d = 1 to the 
exponentially weighted moving average (EWMA). In practice, the difference between the 
Eq. (jsl and the Eq. ( 12 ) is that the function /* places more importance to the data around 
the point ii,i2, id- This is achieved by assigning to the function a weight (1 — a), while 
all the other values are summed together and weighted by a. In Fig. [5} we show the ratio 
'^dma/^^'" obtained by implementing the algorithm respectively with the function / (solid 
lines) and /* (dashed lines) in the range 10 < s < 100. The ratio (t'^ma/ noticeably 
closer to a constant value when the function / is replaced by /*, with a corresponding 
reduction of two orders of magnitude in the relative error AHj. 



V. CONCLUSION 



We have put forward an algorithm to estimate the Hurst exponent of fractals with arbi- 
trary dimension, based on the high-dimensional generalized variance (t]^ma defined by the 
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Eq. Q. 

The methods currently used to estimate the Hurst exponent of high- dimensional fractals 
are based on: (i) 1 — d two-point correlation and structure functions operated along different 
directions, (ii) high— c/ Fourier and wavelet transforms [9l [THl ttH US [131 129] . The advantage 
of the methods [i) is the ease of implementation. Their drawback is the limited accuracy 
due to biases and nonstationarities, being these functions calculated over the entire fractal 
domain. The methods {ii) are more accurate, however their implementation is complicated 
especially for data set with limited extension. The generalized variance ctdma "scaled", 
meaning that it is calculated over sub-arrays of the whole fractal domain by means of the 
function /. The "scales" are set by the size of the sub-arrays rii x n2 x .... x rid. Therefore, 
the proposed method exhibits at the same time: (a) ease of implementation, being based 
on a variance-like approach and (b) high accuracy, being calculated over scaled sub-arrays 
rather than on the whole fractal domain. 

A further important feature of the proposed algorithm is that it can be implemented 
"isotropically" or in "directed" mode to accomplish estimates of H in fractals having pref- 
erential growth direction e.g. biological tissues, epitaxial layers or in crack propagation in 
fracture. The isotropic implementation is obtained by taking 6i = 62 = ... = Od = 1/2 in 
the Eq. Q. This choice implies that the reference point {ii,i2, ...,id) of the moving average 
lies in the center of each sub-array ni x n2 x ... x and thus / is calculated by summing 
the values of / around {ii,i2, Conversely, to implement the algorithm in a preferen- 

tial direction {directed implementation) , the reference point must be coincident with one of 
extremes of the segment ni, or with one of the vertices of the square grid ni x n2 or of the 
d-dimensional array ni x n2 x ... x Ud. The directed implementation can be performed by 
choosing for example 61 = 62 = ... = Od = 0. 

Further generahzations of the proposed method can be envisaged for applications to the 
analysis of time-dependent spatial correlations in > 2. 
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APPENDIX: TRANSFER FUNCTION OF / 



The function /, defined by the Eq. corresponds to the discrete form of the integrah 

1 

f{xi,X2,..., Xd) = / dx[ 

riT2...Td Jxr-Ti 

dx2--- / dx^f{xi,X2,---,x^) (A-l) 

2-T2 -fxd-Td 

where for the sake of simphcity we have considered the case 6i = 62 =,...,= 6^ = 0. 



The Eq. (A.l) can be rewritten as a convolution integrah 

f{xi,X2,...,Xd) = ^ / dxiul — ] I dxiui — 



dx: f/(^^^ fix, - xl, X2 - xl Xd - X*,) (A.2) 



with the convolution kernels given by the boxcar function: 



U{x:/n) 



1 for < X* /Ti < 1 
elsewhere . 
The transfer function can be calculated as follows: 



nriuji,uj2, ... ,ujd) = / dxi dx2... 

nT2--TdJo Jo 

Td 

da;rfexp[-z27r(u;ia;i + UJ2X2 + ... + oJdXd)] (A.3) 



that can be written as: 



1=1 

that is thus (i-times the function smujiTi/ujiTi. 

The Fourier transform JF of the function / can be obtained by means of the following 
relationship: 



T{uji, UJ2, tJd) = Hriuji, UJ2, tud) J'(cJi, UJ2, ...ujd) (A. 5) 

where T{uji,uj2, ..-.ujd) is the Fourier transform of the function /(xi, X2, Xd). 
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The power spectrum S of the function / is given by: 

S{uJi,UJ2,-;UJd) = \Ht{(^1,'^2, ■■■,'jJd)\'^S{uJi,UJ2, ■■■,UJd} (A.6) 

where 8(001,002, is the power spectrum of the function f{xi,X2, ■■■,Xd)- 
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